#include <stdio.h>
#include <math.h>

#include "Hartmann.h"
#include "Util.h"

Hartmann::Hartmann(int n, int m){
	cont = 0;
	this->n = n;
	this->m = m;
}

Hartmann::~Hartmann(){

}

void Hartmann::setFnEvals(int c){
	cont = c;	
}

int Hartmann::getFnEvals(){
	return cont;	
}

bool Hartmann::isNearOptimum(double fBest){
	double bestValue = -1.0;
	switch(n){
		case 3: bestValue = -3.86278;
				break;
		case 6: bestValue = -3.32237;
				break;
	}

	double deltaValue =	fabs(fBest - bestValue);
	double equation;
	
	equation = fabs(bestValue)*0.0001 + 0.000001;	
	if ((deltaValue < equation) || (Util::equals(deltaValue, equation))){
		return true;
	} 

	return false;
}

 double Hartmann::calc(double *x){
	cont++;
	long double value;
  	
	switch (n){
		case 3: value = func43(x);
		case 6: value = func46(x);
	}

	return value;
 }


 // Hartmann function (3,4)
 double Hartmann::func43(double *x){
 	 long double value = 0;
 	 int i, j;
 	 long double qi;

 	 long double A[4][3] = {{3.0,10.0,30.0},{0.1,10.0,35.0},{3.0,10.0,30.0},{0.1,10.0,35.0}};
 	 long double P[4][3] = {{0.6890,0.1170,0.2673},{0.4699,0.4387,0.7470},
				 {0.1091,0.8732,0.5547},{0.0381,0.5743,0.8828}};
 	 long double alpha[4] = {1.0,1.2,3.0,3.2};

 	 for (i=0;i<4;i++){
 	     qi = 0;
 	     for(j=0;j<3;j++)
			qi += A[i][j]*pow(x[j] - P[i][j],2);

    	 value += -1*alpha[i]*exp(-1*qi);
    }

  	return value; 
 }

 // Hartmann function (6,4)
 double Hartmann::func46(double *x){
	long double value = 0;
	int i, j;
  	long double qi;

  	long double B[4][6] = {{10.0,3.0,17.0,3.05,1.7,8.0},
				 {0.05,10.0,17.0,0.1,8.0,14.0},
				 {3.0,3.5,1.7,10.0,17.0,8.0},
				 {17.0,8.0,0.05,10.0,0.1,14.0}};
  	long double Q[4][6] = {{0.1312,0.1696,0.5569,0.0124,0.8283,0.5886},
				 {0.2329,0.4135,0.8307,0.3736,0.1004,0.9991},
				 {0.2348,0.1451,0.3522,0.2883,0.3047,0.6650},
				 {0.4047,0.8828,0.8732,0.5743,0.1091,0.0381}};
  	long double alpha[4] = {1.0,1.2,3.0,3.2};

  	for (i=0;i<4;i++){
      	qi = 0;
      	for(j=0;j<6;j++)
			qi += B[i][j]*pow(x[j] - Q[i][j],2);

      	value += -1*alpha[i]*exp(-1*qi);
    }

  	return value;
 }
